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Abstract 



In the deformation of layered materials such as geological strata, or stacks of paper, me- 
chanical properties compete with the geometry of layering. Smooth, rounded corners lead to 
voids between the layers, while close packing of the layers results in geometrically-induced 
curvature singularities. When voids are penalized by external pressure, the system is forced 
to trade off these competing effects, leading to sometimes striking periodic patterns. 

In this paper we construct a simple model of geometrically nonlinear multi-layered struc- 
tures under axial loading and pressure confinement, with non-interpenetration conditions 
separating the layers. Energy minimizers are characterized as solutions of a set of fourth- 
order nonlinear differential equations with contact-force Lagrange multipliers, or equivalent ly 
of a fourth-order free-boundary problem. We numerically investigate the solutions of this free 
boundary problem, and compare them with the periodic solutions observed experimentally. 

Keywords: bifurcation; localization; void formation; instability; buckling; Kuhn- Tucker 

theory; contact condition. 

1 Introduction 

1.1 Folds in multilayer structures: rocks and paper 

Geological folds arise from the in-plane compression of layers of rocks from tectonic plate move- 
ment, and provide one of the overarching themes of this special issue (see also [ZIII5)- Geometric 
constraints following the need for layers to fit together generate a number of distinctive forms of 
folding, not seen in single- layered systems (see for example Price & Cosgrove [H]). A couple of 
illustrative examples are shown in Fig. [T] On the left is an example of parallel folding where, as 



Figure 1: Folding patterns in thinly bedded shales and sandstones in Cornwall, UK. Left: parallel 
folding at Bude. Right: chevron folding at Millook Haven. 

the eye moves from bottom left to top right, the curvature of the layers increases until eventually 
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it becomes too tight, and forms a distinctive swallowtail singularity [7, after which the layers are 
seen to crumble. On the right, the geometric constraints have been overcome by the rocks adopt- 
ing straight limbs and sharp corners, in what is known as chevron folding [6 . This phenomenon is 
strongly reminiscent of kink banding^ seen also in other contributions to this special issue [131 ES] • 
The sequence of events leading to such structures is of course lost in geological time, but small- 
scale experiments on layered materials such as paper do provide some suggestions. The outcome 
of one such experiment [17 is shown in Fig. [2] About 800 half sheets of A4 cut longitudinally 
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Figure 2: Top: Kink banding in compressed layers of paper. Deformation is forced in about 625 
sheets of paper from a sample of 800. Bottom: output diagram: load vs. end displacement. 



were squashed together under a transverse load between steel plates, and then slowly loaded in the 
longitudinal direction by a second applied axial load. Load cells were used to record both loads 
at one-second intervals, and an in-plane transducer registered the axial end displacement. To aid 
visualization, a black-edged sheet was introduced approximately every 25 layers. 

Output from a typical experiment is shown at the bottom. Kink bands formed in sequence 
from the loaded end, starting at the point B, followed in turn by D, F and G, seen on the output 
diagram. Each band forms with a clear-cut instability, observed by an corresponding instantaneous 
fall in applied axial load and followed by subsequent re-stabilization or lock-up. By comparison. 
Fig. |3] shows various stages in the loading history of a similar experiment, but with soft foam layers 
replacing the upper and lower (stiff) supporting foundation, such that the overburden pressure 
and accompanying stiffness are much reduced. The folds still develop and lock-up sequentially, 
but the sharp corners and straight limbs of the kink-bands are replaced by curved layers, as in the 
parallel folding example of Fig. [T] lock-up is here associated with the tightening of the curvature 
on the inside of a fold [2 . 

The test shown sequentially in Fig. [4] starts under conditions similar to parallel folding, but 
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Figure 3: A compression test with a weaker, but stiffening, foundation (foam). Observe the parallel 
folding and and the serial or sequential buckling. 

evolves to a situation more like kink-banding. Foam is again used as a supporting medium, but 
at the onset of loading only a small overburden pressure is applied, so effectively an Euler buckle 
initially appears over one half wave; this is just apparent in the left-hand image of Fig. [4j The 
long wavelength means that deflections into the foundation grow relatively rapidly, and quickly 
bring the outer steel plates into play, stiffening the foundation. This strongly nonlinear effect both 
flattens the crest of the wave, and instigates two outlying kink bands orientated across the sample 
at about 30° to the transverse direction, seen most clearly in Fig.[4j Although these appear similar 
to those of Fig. |2j they are formed under quite different circumstances, as follows. 

When a kink band is orientated obliquely across the sample, as it displaces from its undeflected 
(flat) state, the layers initially dilate within the band [17 . Those seen in Fig. [5] have used this 
property to select their angles of orientation across the sample, such that maximum available 
dilation exactly cancels the pre-compression induced by the overburden pressure and friction is 
eliminated. They propagate across the sample in this optimum configuration, in which work 
done against friction is minimized; this orientation is then maintained as the band develops and 
eventually locks- up. 

The same is not true for the system of Fig. |4j The development of the band proceeds slowly 
from the onset of loading, and its orientation is not freely chosen but imposed from the start, 
and dictated by factors such as sample length and overall thickness. If the angle of orientation of 
the band is greater than its optimum value, voids will result as seen here. Note that, although 
the introduction of a black-edged sheet clearly helps to induce voids, they appear as a periodic 
sequence every five black sheets, or about once every 125 layers. Very similar effects are also found 
in paperboard creasing and folding, as seen in Fig. 8 of Beex & Peerlings [T (this special issue). 
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Figure 4: A compression test with low initial overburden pressure. Note the periodically positioned 
voids at kink band corners in the final picture. 

1.2 A model of multilayer folding 

In this paper we focus on how bending stiffness, overburden pressure, and geometric constraints 
combine to cause voids, or not. The intriguing periodicity of the voids in Fig. |4] serves as an 
example: is it possible to explain, or even predict, the way voids concentrate at certain layers? 

We do so by formulating a simple model of a multilayer material in which we eliminate all but 
the essential aspects. The layers are described by (linear-elastic) Euler-Bernouilli beams, which 
are inextensible and laterally incompressible; the contact condition is 'hard', and the overburden 
pressure is modeled by a simple global energetic expression (see Section [s] below) . Geometric 
nonlinearity is preserved, and the system will undergo large rigid-body rotations. 

Such a model is highly simplified in comparison to other studies in the literature (see [8 for a 
recent review), and differs from a large part of the literature by assuming elasticity rather than 
viscosity. By studying such a simplified model, however, the appearance and size of voids are 
precisely defined, and we have detailed access to their properties. As for the elasticity, we believe 
that for the type of question that we pose here elastic and viscous materials behave sufficiently 
similarly on the relevant time scales. In other work we have seen how the assumption of elasticity 
allows for a precise analysis [TTl |T7l [2j [3] . 

Our previous study of the behaviour of voids concentrated on a single infinite elastic inextensi- 
ble beam forced into a V-shaped corner under uniform pressure q . We generalize this model to 
multilayer folding in two steps. First we recast the model of voiding of [6] to embrace finite length 
beams and axial loads. The introduction of axial load negates the favourable properties found for 
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minimizers of the purely pressured case: minimizers are no longer necessarily convex, symmetric, 
or have a single void space. It is shown that, under various loading conditions, minimum energy 
can lie with buckling on a local length scale [H [10] , demonstrating higher order accommodation 
structures. 

The main section of this paper then combines the multilayered set up presented in [2 with the 
elasticity properties of Section [2j to capture the periodic behaviour of Fig. [4] A first experiment 
in this direction is to stack identical copies of elastic layers; an elementary geometric argument 
suggests that reduction of the total void space could force such a stack into a straight-limbed, 
sharp-cornered configuration, as seen in Fig. [Sj 




Figure 5: Sharp- angle, straight-limb folds give rise to smaller void space than rounded folds. 



In Section [3] we construct this experiment mathematically by deriving a potential energy func- 
tion for a system of N identical layers that are forced into a V-shaped corner by a uniform axial 



load P and overburden pressure q (see Fig. 11 below), to produce a so called 1-period solution. 
The energy for this specific case maybe simply expressed in terms of a single layer, which leads to a 
similar analysis as presented in Section [2| The results match our initial elementary argument and 
show that the total void space reduces as the overburden pressure increases, giving deformation 
with straight limbs and sharp corners. If we then consider m-period solutions, where one void 
forms every m layers and then repeats, we are led to an interesting question. Which m-periodic 
solution has minimal energy? Constructing an energy for m-periodic solutions and the subsequent 
variational analysis is non-trival, and we therefore leave this to a future publication. We do, 
however, offer some preliminary energy calculations, which suggest that the energy typically is 
minimized at an m that is larger than one but still finite. When we choose reasonable values for 
the parameters in the experiments described in Fig.[4| the predicted value of m is of the same order 
of magnitude as the observed value, giving us the hope there is some more interesting research to 
follow. 



2 Single layer bending and buckling 

Dodwell et al. ^ treats the system of Fig. [g] as an obstacle problem with a free boundary, where the 
potential energy for the system is minimised subject to the pointwise obstacle constraint w > f. 
The potential energy functional, V{w), comprises strain energy stored in bending, Ub^ and work 
done against overburden pressure in voiding, /7y, so that 

V{w) = Ub + Uv^^ J^—'^^dx + q {w-f)dx, (1) 

^ ^-00 V-'- + '^x) J -00 

where B is the elastic bending stiffness of the layer and q is the overburden pressure per unit 
length. Note that Ub and Uy are measured from different zero-energy configurations; the former 
is zero in the flat state and the latter zero in the V-shape. For finite length beams, the two can be 
reconciled by noting that the work done by q in moving from the flat beam configuration to the 
displaced state is a constant (the area of the triangle between the V and the flat beam) minus the 
work done against q in voiding: constants of energy are known not to affect equilibrium equations. 
The same idea extends in the limit to beams of infinite length. 
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Figure 6: Single elastic layer, of infinite length, being forced into a corner by the action of pres- 
sure q. 

By considering stationary solutions of ([T]) with respect to vertical displacement, equilbrium for 
the system is described by the Euler-Lagrange equation. 
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where /i is a Lagrange multiplier, arising from the point- wise obstacle constraint w > f. Physically 
this quantity can be interpreted as the vertical component of the reactive contact load with the 
boundary /; therefore /i is zero wherever there is no contact, it matches the overburden pressure 
q along stretches of contact, and equals the discontinuous shear force at points of delamination. 
Dodwell et al [6 prove that there exists a unique minimising solution to the Euler-Lagrange 
equation and that such solutions are convex and symmetric, and they have a single finite 
interval in which they are not in contact with the boundary. 

Equation ([2| is written above in Cartesian coordinates {x^w), but exists also in intrinsic co- 
ordinates (s,?/^), where s is arc length and ip = dw/ds. Both representations have their uses: 
intrinsic coordinates give simpler equations, but for multi-layered problems they carry the added 
complication of a different coordinate frame for each layer. 

The formulation of Dodwell et al [6] is readily extended to embrace the different configurations 
of Fig. [7| for a long beam of finite length L, with dead load P acting down the slope. Here we 




Figure 7: Finite length beam under pressure loading q and axial compression P. 
are primarily interested in the bottom picture, where the combination of P and q has pushed the 
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beam into a similar configuration as Fig. |6j such that some portion of the length is in contact 
with the supporting medium. The shortening A, compared to the flat state Wx = 0^ and measured 
along the limb, is given by 

L ^ 



A = -£^/rT^+-(v/rT^-i)+ / v'TT^dx, (3) 

^ Jo 

where i is the horizontal non-contact length as seen in Fig. [6j The work done by the compressive 
load P is then PA, and after ignoring a constant the total potential energy over the non-contact 
region can be written, 

V{w) = Ub^Uv-2PA^- / "2W2 ^^ + ^ / i^-f)dx 



(1+^2)5/2 

-2P (^^l^k^ - Vl + wl dx j . 



(4) 



This has the Euler-Lagrange equation. 
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q = 0. (5) 



Note that the effect of the dead load P can be replaced by rigid load more closely resembling the 
experimental set-up, by introducing the extra constraint A = constant. 



2.1 Numerical solutions 

We now consider solutions which minimise V subject to the constraint A = constant. With the 
additional loading constraint the favourable properties of minimisers for the purely pressured case 
(unconstrained minimization) may cease to hold: constrained minimizers of V are not necessarily 
unique, symmetric, or convex. Indeed, although there remains the possibility of higher order, 
asymmetric solutions [121 S] ? we shall here restrict ourselves to symmetric solutions that have a 
single finite interval for which they are not in contact with the obstacle /. 

The Lagrange multipler /i can be eliminated by writing equation ([5| as a free-boundary prob- 
lem [6j. Numerically this problem can then be treated as an initial boundary problem with at 
most two unknown shooting parameters^ i and P, depending on the form of loading. Solutions 
of ([5| are found by shooting from the free boundary at x = — ^ with the initial conditions 

w{-i) = M, Wx{-1) = k, Wxx{-i) = 0, Wx,x{-i) = -(1 + k^f^^^ 

towards the target symmetric section at x = given by 

^x(O) and ^a^a^x(O) 0. (7) 

This initial conditions ([6| can be found by integrating ([5| once, and imposing continuity of w^^ 
and Wxx at points of delamination x = 

Fig. [8] shows a number of numerically obtained solutions of this free boundary problem. 
Fig. JsFa) gives the unique equilibrium solutions for P = and different values of whereas 
Fig. pTb) displays plots of one of the equilibrium solutions for g = 1 and increasing P. 

Fig. [9] shows numerical solutions of ([5|, on a plot of load P against its corresponding deflection, 
as measured by A, at a constant value of q. This shows two of the many possible types of solution. 
From P = 0, we first follow the solid line corresponding to the simple convex solution oflsFb). 



k 



q£ 



, (6) 
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Figure 8: Equilibrium configurations for k = 1. (a) P = and < q < 100. (b) Primary solutions 
for q = l and < P < 50. 




Figure 9: The main graph shows values of P for given A for two types of solutions of the Euler- 
Lagrange equation ([5|, for fixed pressure q = I and angle of limb fx = k = 1. The smaller plots 
show the profile of the layer at the positions indicated. Am marks the position of a Maxwell 
displacement, where energy levels in the two post-buckled equilibrium states are the same [17 : for 
A > Am, minimum energy rests with the upbuckled shape. 



As the layer is forced into the corner singularity we observe that A ^ L/2 and P ^ oo. The 
vertical line in Fig. [9] marks the Maxwell displacement, Am, where energy levels in the two 
equilibrium states are the same p/7j. If A > Am, minimum energy rests with the upbuckled 
shape. In particular, A = L/2 gives a upper bound for the existence of simple convex solutions 
beyond which only higher-order solutions, such as that indicated by the line-dash curve, exist. 
It is worth noting that such higher modes apparently do exist in geological structures, known as 
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accommodation structures [HI p. 381]. Fig. 
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shows one example from a rock outcrop in North 




Figure 10: The formation of higher-order accommodation structures in the hinge regions of a 
multilayered buckle. Rock sample from N. Wales, after Price & Cosgrove [14j. 

Wales. The possibility of energy minimisers involving buckling on a smaller scale is particularly 
interesting, but we leave the investigation of such convoluted solutions to a further contribution [5 . 



3 Multilayered bending and buckling 

3.1 Multilayer formulation 

We now extend the single-layered model to a multilayered formulation. Consider a stack of N 
identical layers of uniform thickness h and length L forced into a V-shaped singularity by an 
overburden pressure gradient defined by q on the top layer, and axial load P per layer, as seen in 
Fig. [Til 

The centreline of the i^^ layer from the bottom is labelled as w'^{x), and the top and bottom 
boundary of that layer as wl{x) and wl{x) respectively. A constant axial load P acts on the stacks, 
which shortens the i^^ layer by a distance parameterized by A^. As before, is measured from 
the centre of the fold, and here additionally gives the horizontal length of the layer. We also 
introduce the obstacle constraint wl{x) < wl'^^{x) for i = 1 to — 1 and w^{x) > /, so that 
adjacent layers cannot interpenetrate or pass through the boundary f = k\x\. 

We now construct a potential energy function for the system, by adding together the contribut- 
ing parts. 

3.2 Bending energy 

The total bending energy of the layer is the sum of bending energy of each layer, given by 

„ f ^ 



The quadratic dependence on w^^ implies that sharp corners have infinite bending energy, and 
therefore for finite overburden pressure, layers will void rather produce singular corners. 
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Figure 11: Multilayered setup. Here N layers are forced into a V-shaped boundary by an external 
horizontal load P and pressure q. 



3.3 Work done against overburden pressure 

In experiments with paper, where voids occur there is clearly a subtle interplay between release 
in overburden pressure in these areas, and some compensatory increase in internal strain energy 
elsewhere. Rather than attempt to model this complex situation directly, we simply assume here 
that the creation of all voids does work against a constant pressure regardless of position in 
the sample. Further development is clearly necessary for application to experiments like that of 
Fig. [4j but we leave this for future work. 

The overburden pressure acting on the layers is q per unit length, and therefore the total work 
done against overburden pressure is q multiplied by the sum of all voids between the layers, given 
by 



r r 

Uy = q {w^-f)dx^qy2 {wl""^ - w\) dx . 



The first term gives the energy associated with the void formed between the bottom layer and 
the boundary /, and the second calculates the size of the voids between each of the remaining 
layers. If q is large, Uy implies a severe energy penalty. 

The fundamental assumption of engineering bending theory, that plane sections remain plane 
and normal to the centreline, means the top and bottom boundaries of a layer can be calculated 
parametrically by propagating the centerline forwards and backwards a distance h/2. In this 
case it is necessary for the bottom of one layer, wl'^^ , and the top of the next, wl, to be defined 
according to the common fixed coordinate system (x^w). At a given x we wish to calculate the 



void between the i^^ and {i + 1)^^ layer, defined as v'^{x) = wl'^^ — w\. Figure Il2| (a) shows 
that w'^^^ can be determined by propagating w^^^ backwards from the point i), and wl{x) by 
propagating w'^ forwards from ii). These centreline positions are not known a priori, but can be 
reverse-engineered from the common position x by the first-order approximation 



-v/TTRF w\^w' + -^l + {wl)\ (8) 



Fig. 12 (b) demonstrates this approximation, where in) and iv) are found by constructing right- 



angle triangles. This is a thin-layer approximation, in which we assume that the layer is much 
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Figure 12: a) Boundary functions determined via front propagation, b) First-order approximation. 



thinner than the radius of curvature of the layer. Whether that is justified or not depends on the 
situation. We also note that the approximation is exact for straight segments, whereas for the 
convex segment shown in Fig. [12] b), Hi) is an underestimate and iv) an overestimate: and vice 
versa for a concave segment. With this approximation the total energy in voiding becomes 



P N-i . , , 



dx. 



3.4 Work done by load 

If we consider a constant axial load P acting on each layer of the multilayered stack, the total 
work done by the load is the sum of the work done by the load in the shortening of each layer. 
Therefore up to a constant the total work done by the load is, 

N N 



Up = 2P^ A, = -P^l + ^'^^i. 

i=i i=i 

3.5 Total potential energy function 

The total potential energy for the complete system can be written as the sum of work done in 
bending by each layer, total work done against overburden pressure in voiding, and work done on 
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the system by the axial load P in shortening, as follows: 
V{wi, ...,Wn) = Ub + Uv + Up 

+ 9 E / -wi)dx-2Y^ PA,. (9) 

i=l '^^i i=l 

Appropriate differentiation leads to a system of N nonlinear fourth order equations, coupled by 
contact conditions. 



4 Periodic solutions 



The formulation given above gives a general energy description for N identical layers forced into 
a V-shaped singularity. Now, motivated by Fig. [4j we restrict ourselves to minimizing V over the 
class of symmetric periodic functions with one void per period. We also make the assumption that 
layers are symmetric, with straight limbs of fixed angle k connected by a convex curved section. 



From Section 2.1 we recognise that higher-order solutions and asymmetric solutions might exist, 
but the experiments shown in Fig. [4] suggest that the assumption of this simplified geometry is 
justified at least in some limited context. 



4.1 1-periodic solutions 

We start with the simplest periodic solution where the structure repeats every layer, hereafter 
named a 1-periodic solution. Since each layer deforms identically the total potential energy may 
be described in terms of a single layer. Each layer does the same work in bending and shortening, 
and the size of the void can be calculated by shifting the bottom boundary of a layer to fit above 
the top of the same layer. Thus, for the primary buckling mode shape of Fig. IsF 



w 



The total potential energy (|9| simplifies to 



V{w) 



BN 



•^2)5/^ 



dx 



+ g J{w-f)dx-q{N-l)h J v/1 + K) 



2 dx - NPA, (10) 



where the shortening A is again given by ([3|. Minimisers of ( 10 ) solve the Euler-Lagrange equation 
over the non- contact region, given by 



BN 



(1 + ^2). 
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q = 0. (11) 



As in Section [2j boundary conditions can be found by integrating ( |TT] ) with respect to the spatial 
coordinate x and imposing continuity of Wx and Wxx at points of delamination x = ±£, so that 



wi-e) = M, wx{-i) = k, w^:,{-e) = 



(12) 



and 



{qh{N - 1) + 2NP) 



k 
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coupled with the symmetric boundary condition at x = 0. Figure [13] shows examples of solution 
profiles obtain by numerical shooting methods. The results match physical intuition in the sense 
that as overburden pressure q or load P increases the void size decreases, leaving in the limit 
deformations with straight limbs and sharp corners. 





Figure 13: Two examples of solution profiles for fixed angle k = 0.75, h 
two values of overburden pressure; (left) q = 10 and (right) q = 100. 



0.02, load P = 1, and 



4.2 m-periodic solutions 

The previous section naturally extends to m-periodic solutions, where a void forms after m layers 
and then the structure repeats. This leads to various interesting questions. Is there a value of m, 
not always 1 or oo, which has minimium energy among the class of such periodic solutions? Can 
we predict the value of m = 125 seen in the periodic solutions presented in Fig [4]? As we have seen, 
1-periodic solutions can be solved using previously developed analysis [6]. This is not the case 
for more general m-periodic solutions, where deriving the energy and the subsequent variational 
analysis is a non-trivial exercise. We therefore leave the details of this work to a future publication. 
To test the basic ideas, however, we perform a very rough calculation for the experiment presented 
in Fig. [4j to see if a simple energy analysis, with reasonable values of the physical constants, gives 
a value of m of the correct order of magnitude. 

The central premise is that in equilibrium the bending energy in a single m-periodic packet is 
approximately equal to the energy required to form one of the voids. For this analysis we assume 



that each layer has straight limbs connected by a circular arc (see Fig. 14) where the circular arc 
of the outside layer is an arc of radius R swept through an angle 20. The angle is fixed for given 
gradient of the limbs k. The bending energy stored in an m-periodic packet of layers is therefore 



Ub = 



EI 



.-1 



{sf dSi = 



EI 



m ^0 

~. J-e 



where Hii is the curvature and Si is the arc length of the i^^ layer. Here we number the layers as 1 
to m, from the outside inwards. Following a result presented in [2 , hCi = /^i/(l — {i — l)hf<ii)^ and 
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Figure 14: Shows the construction of an m-periodic package as straight hmbs connected by circular 
arcs. Here the thickness of each layer is h and the radius of curvature of the out side arc is given 
hyR. 



it therefore follows that 



:=/(/^i) 

For large m, we consider f{i^i) to be a Riemann sum, and approximate it by the corresponding 
integral, i.e. setting Xi = h{i — 1), 



we approximate > by / dx = — log(l — xk,i 



X = Xm-\-l 
X — Xi 



by which 



EI f 

/7b = - ^ J log(l - m/i/^i) dO. 
For a circular arc ki = 1/R is constant over 0, therefore 



Ub = -^IoJi-^)0. (13) 



Estimating the work done to produce a void is simply the overburden pressure q times the size of 
the void, we can equate the forms of energy, so that 

- ^ log (^1 - O^qx Void Area (14) 

The table below displays estimates of the physical constants for the experiment, and gives a brief 
description of their source. 
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Constants 


Values 


Source 


E 


5kNmm~^ 


117 






I 




Standard Calculation 


Transverse Load 


mkN 


Fig. 


I 


1 


Area 


3 X lO'^mm^ 


Area of A 




paper 


q 


1.33Nmm-2 


Transverse Load/ Area 


Void Area 


lOmm^ 


Fig. 


4 







1 


Fig. 


1 




Radius of Curvature 


5 mm 


Fig. 


I 




h 


0.1mm 


Thickness or 




v5 paper 



Subsituting in the parameter from the table above into ( |14[ ) we obtain 

log (l 



mh 



-2, 



giving a rough estimate of m ~ 40. 

The fact that an intermediate value of m is chosen suggsts that the basic premise of the calu- 
lation (balance of bending and voiding energy) might be reasonable. The fact that the predicted 
number m ~ 40 is in the same ballpark as the observed value of 125 in Fig. [4] further reinforces 
this suggestion. 



5 Concluding Remarks 

This paper extends a single-layered model for folding and void formation [6 to the simplest 
periodic multilayered structure, including the vital effects of axial loading. The model extends our 
insight into the interplay between material and loading properties on one hand and the geometric 
constraints of layers fitting together on the other. The results match physical intuition, in that 
as overburden pressure increases, void size decreases, in the limit leading to deformations with 
straight limbs and sharp corners. Although an oversimplification in terms of the paper experiments, 
and even more so of their true geological counterparts, it does capture many of their important 
hallmarks, as well as demonstrating the possibility of periodic solutions arising naturally in an 
elastic context. 

Studying these 1-periodic solutions has led to interesting questions about the possibility of 
m-periodic solutions. In particular, can we find which m-periodic solution has minimum energy? 
The simple energy argument presented in the previous section gives a rough estimate for m for 
the experiment presented in Fig. [4j This gives an interesting path of research to follow, and one 
we wish to consider in a future contribution. 
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